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We present the strongest current constraints on primordial non-Gaussianity (PNG) from a collec- 
tion of galaxy surveys, combining large-scale clustering measurements and their cross-correlations 
with the cosmic microwave background. We update the data sets used by Giannantonio et al. 
(2012), and broaden that analysis to include the full set of two-point correlation functions between 
all surveys. In order to obtain the most reliable constraints on PNG, we advocate the use of the 
cross-correlations between the catalogs as a robust estimator and we perform an extended analysis 
of the possible systematics to reduce their impact on the results. To minimize the impact of stellar 
contamination in our luminous red galaxy (LRG) sample, we use the recent Baryon Oscillations 
Spectroscopic Survey catalog of Ross et al. (2011). We also find evidence for a new systematic in 
the NVSS radio galaxy survey similar to, but smaller than, the known declination-dependent issue; 
by weighting the observed data with randomized catalogs in right ascension and declination, such 
problems can be moderated significantly. We find no evidence of primordial non-Gaussianity; for 
the local-type configuration we obtain for the skewness parameter —37 < /nl < 25 at 95% c.l. 
when using the most conservative part of our data set, improving previous results by 30%; we also 
find no evidence for significant kurtosis, parameterized by Anl- In addition to PNG, we simultane- 
ously constrain dark energy and find that it is required with a form consistent with a cosmological 
constant. 



I. INTRODUCTION 

It is widely assumed that the density perturbations 
that seeded the cosmic microwave background (CMB) 
anisotropies and later induced the formation of cosmic 
structure were produced by primordial quantum fluctu- 
ations in the early universe [Ij. An early inflationary 
phase is then necessary to stretch such fluctuations be- 
yond the horizon, and to address a number of other issues 
related to the initial conditions of the big bang. While 
the general paradigm of inflation is widely accepted, the 
precise details of the physical mechanism that produced 
an inflationary phase in the early universe remain largely 
unconstrained (see Ref. for a recent review). However, 
any speciflc inflationary model is testable, as it leaves ob- 
servable traces in the later evolution. For example, the 
simplest single-field slow-roll inflation model predicts a 
nearly flat universe and near scale-invariance of the pri- 
mordial perturbations, both of which have been verifled 
e.g. by the WMAP satellite ^SnSj. Additionally, the pri- 
mordial perturbations are predicted to be purely adia- 
batic, which is also indicated (see Ref. 57| for recent con- 
straints on isocurvature) . Finally, the distribution of pri- 
mordial perturbations should be nearly Gaussian [51 [5] : 
testing this last prediction is our present focus. 
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Determining the amount of primordial non- 
Gaussianity (PNG) is therefore instrumental in dis- 
tinguishing between models of the early universe TOl. 
This has been typically constrained from the higher- 
order statistics (such as the bispectrum) of the CMB 
[H [11] ; the latest CMB measurements from the nine- year 
analysis of data from the WMAP satellite constrain 
local PNG, as quantified by the skewness parameter 
/nl (defined in Section |ll| to be — 3 < /nl < 77 at 
95% confidence [5l[6|. Recent CMB constraints on the 
kurtosis parameter jinl are —7.4 • 10^ < ^nl < 8.2 • 10^ 
[T^ . Higher-order statistics of the large-scale structure 
(LSS) may also be used [13 [H], but the difficulty of this 
approach is in disentangling PNG from late-time non- 
Gaussianity produced by non-linear structure formation. 
This is also an issue if using the effect on the non-linear 
contributions to the matter power spectrum [15^ . probed 
e.g. by weak gravitational lensing [16|. The abundance 
of massive clusters is also a sensitive, complementary 
probe of PNG [T7H^ . and is less sensitive to later 
non-linear evolution. 

In addition to these established methods, a new tech- 
nique has been discovered and developed in the recent 
years: it has been shown that PNG of the local and or- 
thogonal types alters the biasing of dark matter halos 
and galaxies with respect to the underlying density field, 
making the bias strongly scale-dependent; thus by mea- 
suring the galaxy bias it is possible to constrain PNG 
[^m31j . This has been exploited by some authors [32J- 
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157] , finding constraints which are competitive with CMB 
bispectrum analyses: Ref. |32] found —29 < /nl < 70 
at 95%, while Ref. [35] reported hints of detection at 
8 < /nl < 88 (95%) for the local model. A measure- 
ment of 3nl was obtained from LSS data by Ref. [38] 
assuming /nl = by rescaling the /nl constraints by 
Ref. [31]: -3.5 • 10^ < ^nl < 8.2 • 10^ (95%). The bias 
of galaxy clusters is similarly affected by PNG, and has 
been used in combination with cluster abundances to con- 
strain /nl [39] • The scale-dependent bias technique will 
remain competitive with the arrival of future galaxy and 
cluster surveys, and is potentially more powerful than the 
cosmic- variance limit of the CMB [JDHIS]. 

The large-scale divergence in galaxy bias produced by 
PNG can also be observed through cross-correlations of 
galaxy surveys with the CMB, provided some of the CMB 
anisotropics are created locally. This is conveniently pro- 
vided by the integrated Sachs- Wolfe (ISW) effect [47j . 
where small secondary CMB anisotropics are induced by 
the decay of the gravitational potentials when the Uni- 
verse undergoes a transition from a dark matter to a 
dark energy or curvature dominated phase. While the 
nature of dark energy is still unknown, a broad range of 
observations have confirmed its presence [H], including 
observations of the ISW effect [see Ref. |49l and refer- 
ences therein]. In most realistic models the ISW effect 
is too small to be directly observed in the primary CMB 
temperature power spectrum, but it can be detected by 
cross-correlating the CMB temperature anisotropics with 
tracers of the potentials, such as galaxy catalogs [50] . 
The latest compilations which examine multiple catalogs 
are those by Refs. [FT and ^2] (G08 herein); the latter 
was recently updated in Ref. [49] (G12 herein) and we 
refer to this paper for a more detailed introduction and a 
longer description of the ISW effect and the data. These 
cross-correlations are sensitive to the galaxy bias, and 
therefore also depend on PNG as above. 

In this paper we use an updated compilation of LSS 
and ISW data to make robust and accurate PNG mea- 
surements. We first review the galaxy samples that 
we consider, revising those used by G12. In particu- 
lar, we include the latest available luminous red galaxy 
(LRG) photometric redshift sample [S^, also used for 
ISW measurements by Ref. ^54], which were calibrated 
using data from the Baryon Oscillation Spectroscopic 
Survey (BOSS) CMASS sample. Refs. [H |Ml HT] 
have shown that properly addressing the possibility of 
spurious fluctuations in the observed number of galaxies 
is crucial in order to measure and interpret LSS mea- 
surements. Refs. [371 HSl US] have shown that these is- 
sues are particularly important when measuring PNG. 
We address these concerns in all of our analyses and 
present a new method to reduce the impact of some well- 
known systematic effects present in the radio-galaxies of 
the NVSS survey. 

We then explore the cosmological consequences of 
our data, using a nested-sampling Monte Carlo Markov 
Chain (MCMC) method and marginalizing over extra 



nuisance parameters to account for a range of possible 
uncertainties and systematics. We derive cosmological 
constraints for both PNG and dark energy, extending 
the data set to include both the galaxy-CMB and all 
the galaxy-galaxy correlations, and include their full co- 
variance. Indeed, we demonstrate the importance of the 
cross-correlations between different galaxy catalogs as a 
robust method to estimate the amount of PNG in the 
data. 

The plan of this paper is as follows: we briefly re- 
view the theoretical basis of the LSS in the presence of 
PNG and the ISW effect in Section |ni We then describe 
our data set and the most important systematics in Sec- 
IV] clarifies how we model the ob- 
VTwe discuss the cosmological in- 



tion |III[ while Section 
servations. In Section 
terpretation of our data, before concluding in Section [VT 



II. THEORY 
A. Primordial Non-Gaussianity and the LSS 

We consider here deviations from Gaussianity at the 
three- and four-point level (skewness and kurtosis). For 
PNG of the local type, the Bardeen potential at early 
times (zi,) is assumed to be written in terms of an auxil- 
iary Gaussian potential ip as 

^(x,z*) = </?(x,z^) + /nl [f'^{x,z^) - ((^^)(z*)] 

+gNL ['pHyi,z^)~3{ip^)iz^)ip{x,z^)] , (1) 

where /nl and ^nl quantify the amount of PNG. In this 
work we will also study simple deviations from the local 
case, as described below. 

It was originally shown by Ref. [21] that primordial 
non-Gaussianity induces a strongly scale-dependent bias, 
whose features appear on the largest scales for both 
dark-matter halos and galaxies, and their analysis was 
later clarified and improved [5311571 |311 EOHM]- The 
scale dependence originates through the coupling of long- 
and short- wavelength modes in the density perturbations 
which arise in the presence of PNG by applying the peak- 
background split technique to Eq. ([I]). In Fourier space, 
we can use the linear Poisson equation to show that the 
late-time density can be related to the initial potential 
field by 

S{k, z) c:i a{k, z) Lp{k, Zi,), (2) 
where the tilde denotes Fourier transformation. Here, 
2k^T{k)D{z) .g(0) 



a(fc, z) 



(3) 



where T[k) is the density transfer function, D{z) is the 
linear growth function, g{z) cx (1 -|- z) D{z) is the poten- 
tial growth function, for which g{Zi,) / g{0) ~ 1.306. 
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In the presence of PNG, the Gaussian bias focauss is 
altered by two corrective terms, giving an effective bias 

&off(fc, /nLj^Nl) = &Gauss+^6(/NL,gNL) + A6(A;, /NL,ffNL) • 

(4) 

The first, scale-independent correction (5&(/nl,5nl) is 
produced by the effect of PNG on the halo mass function 
and it is typically negligible compared with the scale- 
dependent part A6(fc, /nl, 5nl); in the following we will 
absorb it into the constant part of the bias, which we 
define as 6i(/nl,5nl) = &Gauss + ^^(/nl, ^nl)- For local- 
type non-Gaussianity, the scale-dependent term is (at 
any given redshift) jMj 



A6'°^(fc,/NL,5NL) 



/?/ /nl + Pg pnl 



(5) 



The coefRcients and /3g will in general depend on red- 
shift and mass (or bias). Here the /nl part is 



(6) 



where 5c = 1.686 is the spherical collapse threshold and 
6i = 6i — 1 is the Lagrangian bias. 

For the (7nl contribution /3g, we use the fit to A^-body 
simulations of the Edgeworth expansion introduced by 
Ref. El: 



(M) -0.7-f 1.4(f>-l)2+0.6 (£^-1)= 

dn'^\M) (v 
dlncr-i 



(7) 



with V = \/Tc 
by Ref. 



.6^ -t- 1, where Sc — 1.42, which was found 
to improve the agreement between the halo 
collapse model and TV-body simulations with respect to 
the Press-Schechter value. We adopt this value in the 
(jNL fitting formulae to be consistent with the derivation 
of Ref. [nil (although we keep the standard value when 
fitting /nl, as the two parts of the model are separable). 
ti'P [M) is the skewness of the density field for the case 
/nl = 1, smoothed at a mass M, and a the r.m.s. of the 
same field at the same scale. For the skewness and its 
derivative we use the fitting functions given in Eqs. (46- 
47) of Ref. [64]. Though their redshift dependence is 
different, both /nl and ^nl cause the bias to become 
scale-dependent with the same power of fc, albeit with 
different amplitudes; this makes distinguishing between 

65 . 



comes with an impor- 
for a narrow range in 



the two parameters challenging 

The fitting formula of Eq. Q 
tant caveat, as it was obtainec 
halo mass and redshift |64j : furthermore, it was shown 
that the effect of gNL changes depending to the halo oc- 
cupation distribution of the used tracers, not only on the 
Gaussian bias, and this effect becomes more severe for 
less biased samples, such as in our case. This formula 
is, however, the best currently available, and the form 
of the constraints obtained using it will be indicative of 
those that would be available using a more accurate fit. 



We therefore adopt the formula, testing for specific non- 
Gaussian deviations of this form, making sure that when 
inferences are drawn from the results we remember these 
might be different from true gNL-hke deviations. 

Scale-dependent bias was first used to constrain /nl 
from the LSS by Refs. [32l |6T] with varying results, pos- 
sibly caused by different assumptions for the bias evo- 
lution with redshift. Ref. [38_ obtained bounds on (^nl 
from the same data consistent with Gaussian initial con- 
ditions. Some later work [34lf36] found hints of detection 
of a positive /nl, while more recent analyses that take 
into account possible systematic errors in the large-scale 
clustering measurements have found measurements con- 
sistent with no PNG [37] . 

Similar expressions have been derived for other vari- 
eties of PNG, such as orthogonal and equilateral [5^163] . 
featuring weaker or no scale dependence. Furthermore, 
some models also predict a scale dependence of /nl 
[661 167] , often parameterized with a spectral index n f^-^ : 
recently, the first constraint on this quantity has been 



found 



= O.SOtH at 95% c.l. from CMB data IMI. 



Given the difficulty in distinguishing between different 
models with current bias measurements, and to be as 
general as possible, one can adopt a simpler phenomeno- 
logical approach by allowing the scale-dependent term to 
have an arbitrary amplitude and scaling index. To keep 
the interpretation of the parameter constraints closer to 
the usual local case, we adopt a parameterization 



A6s°"(fc,/, 



NL) 



Q:'(fc,aNL) ' 



(8) 



where we have only one extra parameter onl, changing 
the slope of the bias correction: 

, 2k''''^T(k) D(z) 5(0) 

The local, scale-independent case is recovered for onl = 
2. Notice also that this can be related to the local spectral 
index as 

fc""-=fcM^j , (10) 

where fcpiv is the pivot scale. We will use fcpiv — 0.02 /i 
Mpc~^ in the following as in Ref. (42] . 



B. The ISW Effect and Large-Scale Structure 

For a recent, detailed summary of the ISW theory and 
its observational status, we refer to Ref. [JS]. Briefly, if 
the gravitational potentials ^ are evolving in time, 
they will imprint secondary temperature anisotropics 
Q = 6T/T in the CMB in a direction n given by 

e(n) = - y e-^^") ($ + [77, n(7,o - v)] dv , (H) 
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where 77 is conformal time, the dots are derivatives with 
respect to it, and r is the optical depth of CMB photons. 
Similarly, for a galaxy catalog i with redshift distribution 

(pi{z) = dNi/dz{z), its expected overdensity is 

(5g,(n)= J bg,{z)tpi{z)S{n,z)dz , (12) 

where bg. (z) is the galactic bias and 6{n, z) the dark mat- 
ter overdensity. 

From these relations, the auto- and cross-spectra in 
harmonic space can be constructed assuming linear the- 
ory as 

Cf^^ ^'^ j dkk' P{k)Wj>^{k)Wf'{k) 

Cj^s. ^ ^ J dkk'^ P{k) Wl{k) Wf' {k) ; (13) 

here P{k) is the linear matter power spectrum at z = 0, 
and the sources for galaxies and the ISW temperature 
anisotropics are (assuming $ = vj/) 

Wf'(fc)= j dzbg^{k,z)ip,{z)D{z)ji[kx{z)], (14) 

where x is the comoving distance and ji the spherical 
Bessel function. These spectra are related to the auto- 
and cross-correlation functions (ACF and CCF), defined 
as 

K;^«^(^)^(<5,.(n)e(n')), (15) 

where the averages are carried over all pairs at a distance 
d = |n — n'|. They are derived from the spectra by 
Legendre transformation: 

00 

w{d)^ —Y,{2l^\)Ci Pi [cos(z9)] , (16) 
where Pi are the Legendre polynomials. 

III. DATA AND SYSTEMATICS 

It is well known [37l [58j [59] that the large-scale clus- 
tering of the LSS is easily contaminated by systematic 
distortions in the galaxy samples caused by observational 
issues, which may significantly bias the PNG results. In- 
deed, any spurious long-wavelength contaminant, pro- 
duced in the Galaxy (extinction-causing dust, stars, syn- 
chrotron emission) or even in the atmosphere (seeing, 
airmass) will generally introduce extra large-scale power, 
which may naively be attributed to PNG. This issue only 
has a limited impact when studying the ISW correlation 
with the CMB, as there is no particular reason why such 



systematics should correlate with the primordial temper- 
ature fluctuations; however it has the potential to sub- 
stantially bias any result using the large-scale angular 
auto-correlation of these galaxies. However, as in the 
case of the galaxy-CMB correlations, cross-correlations 
between independent galaxy catalogs will be in general 
less subject to such systematics, and can therefore pro- 
vide estimates of the large-angle clustering that are more 
robust than the auto-correlations. 

To minimize these effects, and to obtain robust mea- 
surements of PNG, we study a number of data sets previ- 
ously considered in G12, including the main galaxy, LRG 
and quasar samples from the Sloan Digital Sky Survey 
(SDSS) 169;, as well as independent infra-red, radio and 
X-ray surveys. We examine a number of these in greater 
depth in order to address the impact of potential observa- 
tional systematics on the large-scale clustering measure- 
ments. 



A. Data and General Considerations 

Below we discuss potential contaminants to our LRG, 
radio and quasar samples. The other data sets are un- 
changed from G12: The 2MASS infra-red catalog [7D] 
contains 415,459 galaxies after masking, has a median 
redshift z = 0.09 and a bias b — 1.3 if it is assumed con- 
stant. The final imaging SDSS Data Release 8 (DR8) 
main galaxy sample |71j contains 30 million galaxies at 
z — 0.3 and b = 1.2, while the HEAO survey of the 
X-ray background [75] has a broad redshift distribution 
of estimated z = 0.9 and b = 1.0. See G12 for further 
details. 

We pixellate all maps using the Healpix [73] scheme 
at low resolution A'sidc — 64, corresponding to a pixel size 
of ^ 50 arcmin, since we are primarily interested in large 
scales. We measure the angular two-point correlation 
functions between two pixellated maps a, b of weighted 
masks f"",/^ using the estimator 

(11-0(1-') <'^> 

where — ^ ■ ■ /""fj is the weighted number of pixels 
at separation i?, n indicates the mean expectation value 
in each pixel and A^pix is the total number of pixels in a 
map. This generalizes Eq. (10) of G12. 

Note that the overlap between catalogs is significant, 
both in redshift and in footprint: this leads to significant 
cross-correlations and covariances in most cases, as also 
discussed in G08, G12. 



B. The BOSS CMASS Galaxies 

The MegaZ LRGs [71H7B] used in G08 and G12 have 
an excess of power on large scales with respect to the 
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ACDM model [75]. It was shown by Ref. [S3] that the 
Hkely cause for this excess is a systematic contribution 
to the clustering from stellar contamination. Here, we 
have decided to instead use an updated catalog of LRGs 
at redshifts z ~ 0.55, constructed from the final imaging 
DR8 of the SDSS [53], where the systematics are better 
understood. We use the DR8 photometric redshift (pho- 
toz) sample rather than the data release nine (DR9) [77] 
spectroscopic redshift sample recently used for multiple 
cosmological analyses [371 178H84] , because of the larger 
cosmological volume covered. Even though it relies on 
photometric rather than spectroscopic redshifts, it cov- 
ers a usable area of ~9,000 deg^, rather than 3,300 deg^ 
for the spectroscopic sample, giving an increase in volume 
that allows better large-scale clustering measurements. 



1. The Catalog and its Correlations 

The DR8 LRG sample was selected using the same 
color and magnitude cuts as the SDSS-III BOSS 
'CMASS' (constant mass) sample, defined in Ref. (55] . 
using imaging data from the entire DR8 footprint. Here 
we will denote it as the CMASS sample. It used more 
than 100,000 BOSS spectroscopic redshifts for training to 
estimate photometric redshifts and evaluating the prob- 
ability that an object is a galaxy. (3% of objects that 
are selected are stars.) A detailed analysis of possible 
systematic uncertainties was performed by Ref. [53], in- 
cluding seeing effects, sky brightness and possible stel- 
lar contamination. The most serious systematic problem 
found with present photometric redshift catalogs from 
SDSS imaging was the loss of faint galaxies around stars 
(even relatively faint stars to r ~ 20). Ref. [S3] cor- 
rected this issue, which clearly illustrates the need to be 
diligent about stellar contamination, particularly when 
cross-correlating with other data sets which may also in- 
clude some contamination from galactic sources. 

We use CMASS LRGs with 0.45 < ^phot < 0.65 and 
color term (see Eq. 13 in Ref. [53]), C||> 1.6. This is an 
identical sample to that used in Refs. [SB] US] [SS] , except 
that we do not further divide the sample by photomet- 
ric redshift. Our map exploits the "Agtar" method (see 
Section 4.1 of Ref. [53]) to correct for the systematic ef- 
fect of stars. We also correct for the offset between SDSS 
photometry in the North and South Galactic caps [87] 
(NGC/SGC from here on) using the method applied to 
obtain the "ASouth" results of Ref. [53] . We use the same 
mask for this sample as constructed in Ref. [56 , but use 
a marginally more conservative Galactic extinction cut 
of Ar < 0.18 (compared to the fiducial < 0.20), for 
consistency with our other data sets and previous analy- 
ses. This yields a sample with 801,226 LRGs occupying 
a footprint of 8,891 deg^. The redshift distribution of the 
sample is centered around z ~ 0.5, as shown in Fig. [T] 
By fitting the ACF, we have confirmed the bias b ~ 2.1, 
found by earlier work [53] . 

We calculate the auto- and cross-correlation functions 
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FIG. 1. Redshift distributions of the different LRG samples. 
The CMASS DR8 data are less peaked, but the difference in 
the distribution is small for the ISW purposes. 



of these LRGs and the WMAP7 internal linear combina- 



tion (ILC) map using the estimator of Eq. (17), binning 



the results in the range < -d < 12 degrees. We also es- 
timate the full covariance matrix using the Monte Carlo 
method MC2 of G08. We see the results in Fig. [2] where 
we also compare them with the correlations of the MegaZ 
data. While the ACFs are comparable at small scales, 
the excess large-scale power is significantly reduced in 
the CMASS data. Furthermore, the correlation of the 
CMASS data with the CMB is much lower than for the 
MegaZ, and on most scales comparable to the theoreti- 
cal prediction for a WMAP7 best-fit ACDM cosmology 
(assuming a bias b = 2.1, and no PNG), remaining well 
within the (highly correlated) l-cr error bars. Any signif- 
icant deviations could reveal either residual systematics, 
or PNG. 



2. Discussion 

The differences between the CMASS and MegaZ CCFs 
seen in Fig. [2] are large compared to their errors, and it is 
interesting to investigate the reasons for this discrepancy 
in more detail. 

North-South Divide The footprint of the CMASS 
data is larger than that of the MegaZ data (8,891 deg^ 
compared to 7,400 deg^ after masking), with most of 
the additional area resulting from imaging data in the 
SGC [7T]. As a first check, we tested whether the differ- 
ence in the CCF signal is produced by the different sky 
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FIG. 2. Auto-correlation functions and cross-correlations 
with WMAP (top and bottom panels respectively) for the 
LRG data sets. The 1-a error bars are derived from Monte 
Carlos and are highly correlated; the data points have been 
slightly displaced in i9 for clarity. Dark and light blue data 
show the CMASS and MegaZ LRGs, and the relative theoret- 
ical expectations from the WMAP7 best-fit cosmology (with- 
out PNG) and constant bias of 2.1 and 1.7 are drawn in dark 
and light green (ACFs negative at > 5 deg). The red lines 
show the measured correlation functions for the GMASS data 
when the DR7 footprint mask is applied, while for the ma- 
genta line we have also removed the systematic correction for 
the effect of faint stars [53] . A part of the differences between 
the LRG samples can be accounted for by these two effects. 



coverage of the data. For this purpose, we nieasured the 
CCF of the CMASS data in the DR7 region only. We 
found that the signal increases significantly, as shown in 
Fig. [2] (red line), although it still does not match the 
MegaZ result, which used DR7 (i.e., mostly the NGC 
data). 
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FIG. 3. North-South divide in the CMASS cross-correlation 
data. The results for the South galactic cap are contrary 
to the expected ISW signal, but the cosmic variance for this 
small area is large and it is potentially more subject to ex- 
tinction systematics. 



We next examined whether the NGC and SGC foot- 
prints of DR8 give compatible ISW signals. As we can 
see in Fig. [3j this may not be the case, as the NGC of 
the DR8 data gives results which are in better agreement 
with the MegaZ catalog, but the SGC shows little or neg- 
ative cross-correlation signal, and thus when added to the 
NGC data, lowers the overall CCF. However, given the 
relatively small coverage of DR8 in the SGC (1,300 deg^), 
we should be wary of over-interpreting the results as the 
cosmic variance error on the ISW measurement is large. 
The difference between the NGC and the full DR8 data is 
more significant, but well within the expected statistical 
fluctuations given the size of the highly-correlated error 
bars in Fig. [3j 

Stellar Correction We have also checked whether 
the new correction for the systematic effect of faint stars 
introduced by Ref. [53] can better explain the differences 
between MegaZ and CMASS. We can see from the ma- 
genta line in Fig. [5] that indeed, if this correction is not 
applied to the CMASS data, the resulting CCF becomes 
even higher, close to the MegaZ result. However, even 
this effect is not enough to fully reproduce the earlier 
result on all angular scales. 

Frequency Dependence An important property of 
the ISW signal is that it is expected to be independent of 
CMB frequency, making any frequency dependence an in- 
dication of possible systematics. To quantify the level of 
contamination of the two LRG catalogs, we investigated 
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FIG. 4. Frequency dependence of the CCFs, for the MegaZ 
and CMASS data. We can see that the latter data are more 
stable across the whole frequency range. 



in detail the frequency dependences of their CCFs, as 
shown in Fig. |4j Here we can see that the MegaZ result 
is more frequency dependent, while the CCF from the 
CMASS dataset is remarkably constant for all WMAP 
frequency bands, including the K and Ka WMAP bands 
which are most affected by residual galactic emission. 
This robustness, unmatched in any of the other galaxy 
catalogs from G12, is fully consistent with the ISW in- 
terpretation of the observed cross-correlation signal and 
appears to confirm that the stellar contamination in the 
CMASS data is negligible. From this test, we conclude 
that the CMASS catalog is the most robust data set avail- 
able for our purposes. 

Extinction As discussed in G08 and G12, dust ex- 
tinction can be a major source of systematic uncertainty 
in ISW measurements. This is due to reddening from 
our Galaxy altering the inferred large-scale distribution 
of galaxies, which then introduces spurious correlations 
with the CMB, which may also have residual galactic 
emission. For this reason, we have carefully checked the 
effect of extinction on the measured CCF. 

Extinction can be corrected by excluding the most af- 
fected regions from the analysis. In Ref. l53|, areas with 
a reddening Ar > 0.20 magnitudes were excluded, while 
the cut in G12 was marginally stricter, > 0.18 and 
we have applied the latter in the present analysis. Fig. [5] 
shows how the CCFs vary when lowering the amount of 
reddening which is tolerated. When this cut is made 
stricter, the measured signal increases both for the full 
DR8 and NGC-only areas. The average extinction is 




< 0.08 theory CMASS 



-+- 



-+- 



-+- 



South only 




4 6 8 

iS [deg] 

FIG. 5. Extinction dependence of the CCFs for the CMASS 
LRGs. The top panel shows the result for the whole surveyed 
area, while the middle and bottom panels present results from 
the NGC and SGC only respectively. The black lines refers 
to our baseline cut at Ar < 0.18. The greater variation in the 
SGC area reflects both its smaller area and the higher average 
reddening in this region. 



higher in the SGC; the average in the whole footprint, 
calculated before applying the extinction cut itself, is 
Ar = 0.12; this reduces to A^ = 0.10 in the NGC and 
increases to Ar = 0.18 in the SGC. 



3. Conclusion 

From the analyses performed in this section, we con- 
clude that the CMASS data set is the most robust avail- 
able for large-scale structure and cross-correlation mea- 
surements, as it appears to have low stellar contamina- 
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tion [53]. We also observe that the lower CCF from 
the CMASS data is partly due to a lower CCF in the 
SGC portion of the DR8 data; while consistent with cos- 
mic variance, this could also be potentially connected to 
higher extinction effects in the SGC region. 



C. NVSS Systematics 

The NRAO VLA Sky Survey (NVSS) fSSl produced a 
catalog of radio-galaxies covering the full sky at declina- 
tions S > —40 deg. The clustering properties of these 
sources have been measured by several authors [89ti92j . 
who found excess power on the largest scales compared 
with the ACDM predictions and the analyses of other 
data sets. Its combination of large sky coverage and red- 
shift depth makes the NVSS well-suited for the measure- 
ment of ISW cross-correlations and this has been one of 
the most studied data sets in this context [55Hl00j . Fur- 
thermore, these data have been used for other studies, 
such as to detect CMB lensing [101] and to constrain 
inhomogeneous models |102j . 

Unfortunately this catalog is affected by some well- 
known issues. There is a large uncertainty in the redshift 
distribution of the sources, which can make its cosmolog- 
ical interpretation ambiguous. For many applications, 
this is not a serious issue; for example, it was shown in 
G12 that different assumptions for the distributions do 
not significantly affect the ISW measurements, given the 
size of the expected errors. 

More worryingly, the sample contains clear variations 
in the density of sources with declination that are sig- 
nificantly greater than those expected for a statistically 
isotropic universe [521 US] ■ In particular, fewer sources 
are detected at 6 < —15 deg, as can be seen in the upper 
panel of Fig. [6] In addition, we have found a similar al- 
beit smaller issue in right ascension: the number density 
has a roughly linear decrease, reaching a minimum at r.a. 
> 250 deg; this can be seen in the lower panel of Fig. [6j 

The main effect of such systematics, and of any long- 
wavelength noise in the data, will be to add spurious 
signal on the large-scale power spectrum (or ACF). In- 
deed the presence of spurious large-scale clustering power 
above an assumed cosmological model is often taken to be 
a clear indicator of uncorrected systematics in the data; 
however, we must be careful not to do so here as we are 
trying to measure, or limit, a large-scale PNG signal from 
the data. 

In previous analyses, a number of methods have been 
used to correct the declination systematics: in some real- 
space analyses [531 1111102, including G08 and G12, the 
data were subdivided in declination bands, and expected 
densities were rescaled to match each. Other authors 
[551 [HI] found that this problem was minimized by 
imposing a strict flux cut of F > 10 mJy; this however 
discards a large fraction of sources. Finally, others jlOlj 
working in harmonic space addressed the issue by as- 
signing infinite variance to the m = modes of the data. 
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FIG. 6. Systematic density fluctuations in NVSS. The upper 
panel shows the strong change as a function of declination, 
while the lower panel presents a weaker deviation we have 
found as a function of right ascension. The data have been 
smoothed with a Gaussian beam of r.m.s. cr = 5 deg. The 
colored band are the Poisson 1 — cr regions. 



effectively discarding any azimuthally-invariant mode. 

The method of fitting sub-bands in declination effec- 
tively assumes that any declination offsets are smooth 
and fit by the particular distribution of bands chosen. 
Going beyond the earlier analyses in GOB and G12, we 
have found that the signal fluctuates depending on the 
size chosen for these bands, leading to an uncertain auto- 
correlation. This suggests that we need a more robust 
technique for correcting for these effects if we are to 
use this sample to robustly measure PNG. However, the 
cross-correlation with WMAP is fairly stable to such cor- 
rections: we show below in Section W Al that the G12 re- 
sults on the significance of ISW detection do not change 
significantly. 

We introduce a new method to deal with these system- 
atics, which is inspired by the way in which the BOSS 
collaboration created random catalogs that were used in 
order to determine their clustering statistics. Lacking an 
accurate a priori model for the redshift distribution of 
the survey, random catalogs were constructed by taking 
random points within the angular mask of the survey 
and assigning redshifts randomly drawn from the true 
observed redshifts [57]; such maps were shown to pro- 
vide accurate means of reconstructing the correct large- 
scale power with little bias. For NVSS, it is the angular 
mask, rather than the radial mask, which is unknown. 
We observe fluctuations in the number density that de- 
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pend strongly on r.a. and declination S which are clearly 
spurious; we are unable to isolate the origin of these fluc- 
tuations, so we do not know, a priori, the expected mean 
distribution of NVSS galaxies on the sky. We will assume 
the fluctuations we observe are associated with flux cal- 
ibration errors and that they are separable functions of 
r.a. and 5. (Allowing an arbitrary distribution would 
negate any structure measurement.) 

First we present a procedure to limit the impact of fluc- 
tuations that are only dependent on declination. In this 
case we estimate the mean angular mask by generating 
a random map by assigning a number of random r.a.'s 
(within the mask) to each observed 6. Pixelizing, we ob- 
tain a map, Ri, which has the same average declination 
dependence as the observed map. This can be normal- 
ized to give the expected mean number density for every 
pixel. We then obtain the corrected clustering replacing 
the standard estimator of Eq. (17) with the formula: 



1 

1% 



E 



ud R 



1) hfj^ 

(18) 

where nx indicates the mean density of the catalog X. 
This is equivalent to assuming that the probability of ob- 
serving an NVSS source at a given r.a. modulates exactly 
as plotted in the lower panel of Fig. [6j 

To extend this to the case where we have a systematic 
in both directions, we produce a second random catalog 
R2 constructed with the observed r.a. distribution and 
random declinations. The average number density is as- 
sumed to be proportional to the product of these two 
random maps, yielding an effective estimator 
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no Rl,iR2 
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(19) 



In practice, the random catalogs use ten times as many 
points as are observed, which appears sufficient to remove 
spurious shot noise in the estimate of the mean. In par- 
allel with Ref. in?]: we call this procedure 'shuffling' as it 
shuffles galaxy positions when creating the random cat- 
alog providing the expected galaxy density at any point. 
The extra terms in Eq. ( 19 ) can be considered as empir- 
ical corrections in the r.a. and declination directions to 
the expected galaxy density. 

The effect of the 'shuffling' procedure described above 
on the NVSS ACF is displayed in Fig. [7| We see that 
the raw NVSS ACF, when no corrections are applied, is 
large (red crosses) . This result is comparable to what we 
see when the spatially-dependent systematics of Fig. [6] 
are added to the mocks (orange dot-dashed line in the 
bottom panel). 

The shuffling method works to reduce the power in 
both real data and mocks with added systematics: in the 
former case (top panel), shuffling in the declination di- 
rection only (green triangles) yields an ACF comparable 
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FIG. 7. Large-scale clustering of NVSS and systematics sub- 
traction. Top panel: We plot with red crosses the ACF of 
NVSS without any correction, which presents a large excess 
power on these scales. We then show with blue squares the 
measurement using a declination striping correction as ap- 
plied by Refs. |93l I94| . This largely agrees with the result 
obtained by shuffling in declination only using the estimator 
of Eq. ( 18 1 (green triangles). We finally show with dark blue 



circles the measurement obtained applying the estimator of 
Eq. (19 1, i.e. correcting for systematics in both coordinates, 
which is closer to the ACDM theoretical predictions with no 
PNG (purple, solid). Bottom panel: Here we validate the 
shuffling method with mocks. The data are replotted from 
the top panel, while the lines show the mean of Gaussian 
mocks (with added shot noise) under different assumptions. 
We show in magenta (long-dashed) the case where the mocks 
have been given systematic density fluctuations of Fig. [6] then 
subtracted by shuffling. For the brown (short-dashed) curve, 
shuffling is applied but the systematics are not introduced. 
The orange (dot-dashed) case has the full systematics and 
no shuffling, while the green (triple dot-dashed) case has the 
systematics, and only shuffling in the declination direction. 



with what was obtained with the striping method (blue 
squares) — a similar ACF is obtained with the mocks 
with systematics (green triple-dot-dashed line in the bot- 
tom panel). When the full shuffling in both directions is 
applied to the data, the result is lowered further (solid 
navy circles); a similar effect is again seen in the mocks 
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FIG. 8. Large-scale clustering of the quasars and stellar con- 
tamination. We plot in red (solid) the prediction for a model 
with no stellar contamination and no PNG. The green (dot- 
dashed) and blue (dashed) lines correspond to models with 
stellar contamination and with PNG, which are both able to 
explain the observed excess power (black data points). Con- 
stant linear bias is assumed. 



FIG. 9. Top panel: The angular auto-correlation function, 
ACF, multiplied by the angular scale i?, when the DR6 quasar 
sample is sub-divided based on i-band magnitudes. The solid 
black line is the prediction for the standard ACDM model. 
Bottom panel: The measured ACF when the DR6 quasar 
sample is split based on photometric redshift, z. 



with systematics (long-dashed magenta line) in the bot- 
tom panel. Incidentally, this result is similar to what 
is obtained if shuffling mocks which do not include the 
systematics (brown, short-dashed). 

Finally, we expect that shuffling should remove some 
clustering strength even if there are no systematics 
present, as any true fluctuations in number density that 
align with the r.a. and 5 directions will be smoothed. 
We quantify the expected amount by applying shuffling 
to our mock catalogs, and find that a constant suppres- 
sion factor q ~ 2/3 is a reasonable approximation. The 
fiducial theoretical ACF is displayed using a solid purple 
line, while the mean ACF of the mocks when shuffling 
is applied is shown with a magenta dashed line (bottom 
panel). A factor 1/q has already been applied to the 
mocks to account for the power removed by shuffling: we 
can see there is substantial agreement on these scales. 

We apply the shuffling correction to the NVSS mocks 
when estimating the covariance matrix, also including 
the observed systematic fluctuations in the mocks. As 
described above, shuffling reduces the amplitude of the 
ACF by a factor g ~ 2/3: we account for this bias when 
modeling the NVSS ACF by simply multiplying the fidu- 
cial theoretical ACF by 1/q (and any NVSS CCF by 

1/^/9)- 



D. Quasar Systematics 

Quasar samples are promising for future PNG mea- 
surements, as they are easily observed to high redshifts 
and are highly biased tracers of the matter distribution. 
The DR6 photometric quasar catalog has been used pre- 
viously to constrain PNG, and Ref. [33] found that these 
data alone yield /nl — 62 ±26 (Icr). However, this result 
is in conflict with Ref. [31], who found /nl = (for 
their fiducial 'QSO' case): they used an earher sub-set 
of the DR6 data from the SDSS DR3 [lUg. Ref. ES] 
thoroughly tested the DR6 quasar sample against poten- 
tial contaminants, and found significant effects associated 
with stellar density, the stellar color locus, airmass, and 
seeing. They concluded that the sample was not fit to 
use for /nl measurements. Here we further investigate 
this issue, and argue that the quasar data should only be 
included through cross-correlations with other surveys. 
Future spectroscopic quasar samples will be immune to 
many of these issues. 

The fiducial photometric quasar sample we use is 
drawn from the SDSS data release six (DR6) tM]. The 
quasar catalog |105j contains more than one million ob- 
jects selected via a non-parametric Bayes classifier - ker- 
nel density estimate (NBC-KDE) method, including both 
UV-excess (UVX) and high-redshift samples, up to a 
magnitude limit of i = 21.3. The total footprint is 8417 
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FIG. 10. Observed relationship between the number density of DR6 quasars with respect to the average (riq) and a choice 
of systematics: airmass, sky brightness, galactic extinction in the r-band, and stellar density, in different redshift ranges. 



deg^. Nevertheless, to reduce contamination as much as 
possible, we restrict the analysis to the UVX sample, as in 
G08 and G12, and we reduce the footprint with a stricter 
extinction cut, discarding areas with Ar > 0.14, so that 
we actually use 502,565 sources covering 6,912 deg^. 

Photometric quasars are particularly prone to stel- 
lar contamination, as their faint and compact nature 
matches those of stars. If we assume a stellar contami- 
nation fraction given by k, the expected auto-correlation 
function will be given by 

^obs = ( 1 _ 2 ^qso (^) + ^2 u;-*^-- (^) , (20) 

where we estimate the stellar ACF w^^'^'^{d) by measur- 
ing the clustering of stars from the SDSS survey. The 
distribution of stars over the DR6 footprint is distinctly 
anisotropic since they trace the structure of the Galaxy, 
and the auto-correlation of stars is therefore significant 
at large scales. A residual contamination of the order of 
~ 2.5% could explain the observed plateau in the ACF 
as shown in Fig. [8] as could /nl ^ 50. We found from 
both the local relationship between quasar density and 
stellar density and also by directly measuring the cross- 
correlation function between the quasars and a catalog 
of SDSS stars that the actual contamination fraction is 
only ~ 1% in these data. The quasar ACF data thus 
suggest there are either further spurious fluctuations in 
the quasar density field, or /nl is non-zero (as suggested 
by Ref. [36 ). 

It is reasonable to assume that sub-samples at fainter 
magnitude should be less reliable and more severely af- 
fected by systematics. We therefore measured the ACF 
with cuts on z-band magnitude and we display the results 
in the top panel of Fig.[9| Indeed, the ACF of the faintest 
data (z > 20.5; black triangles) displays much larger clus- 
tering amplitudes than the brighter sub-samples (and the 
difference is inconsistent with the change in linear bias 
and a non-zero /nl)- However, we are unable to isolate a 
particular magnitude at which we find that the system- 
atic relationships become negligible and we do not find 



that the /nl constraints from the DR6 quasar ACF be- 
come stable for samples brighter than any given i-band 
magnitude. 

The bottom panel of Fig. [9] displays the ACF when we 
split the quasar sample by photometric redshift. Both 
the z > 2.2 (black triangles) and z < 1.0 (green stars) 
samples display significant large-scale clustering, that we 
can attribute to, in large part, significant stellar contam- 
ination in these subsamples (confirmed both via cross- 
correlation and local relationships with stellar density). 
However, we find that cutting by photometric redshift 
can introduce systematic fluctuations, as would be the 
case if there were any relationship between the estimated 
redshift and a potential systematic. This is illustrated by 
the fact that the ACF of both the 1.6 < z < 2.2 (blue cir- 
cles) and 1.0 < z < 1.6 (red squares) samples has larger 
amplitudes than the z < 2.2 (cyan triangles) ACF, de- 
spite the fact that the z < 2.2 sample contains quasars 
with z < 1.0 and thus has significant stellar contamina- 
tion. Thus, while the 1.0 < z < 2.2 sample removes data 
with stellar contamination, the cuts on redshift appear to 
induce other systematic fluctuations. Some examples are 
shown in Fig. |10[ e.g. we find a significant relationship 
with Galactic extinction Ar, when we cut the sample to 
1.0 < z < 2.2, but find no such relationship when no 
redshift cut is applied. We further find significant corre- 
lations with air mass, sky brightness, and especially at 
high redshift, with the surface density of stars. 

We have not studied the DR6 quasar sample as thor- 
oughly as Ref. [59], but we have, independently, reached 
a similar conclusion: the DR6 quasar ACF should not 
be used to obtain /nl constraints. We have found that 
significant correlations exist with potential systematics, 
such as stellar density. Galactic extinction, and airmass, 
exist in excess of those expected. Also, these relation- 
ships depend non-trivially on the redshift and magnitude 
of the quasar sample that is selected. The effect of /nl on 
the quasar density field will also depend on redshift and 
the bias of the sample, and we are therefore unable to 
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FIG. 11. Complete set of the two-point functions we use. The top row shows the CMB-galaxy correlation functions, while the 
remaining panels are the galaxy-galaxy correlations. Error bars are from 10,000 Monte Carlos, whose means are the red dashed 
lines, and the blue line is the standard ACDM cosmology from WMAP7, with constant biases (not a fit to these data). 



address systematic concerns using the methods outlined 
in Refs. [53, 56 . 

However, we do not expect these issues to be corre- 
lated with other samples, and should be able to trust 
correlations between the quasars and other data sets. In 
particular, the quasars have a large overlap in redshift 
with the NVSS data. Potential SDSS systematics, such 
as airmass and seeing, are survey-specific and should thus 
have no correlation with NVSS data. In addition, we find 
no correlation with NVSS data and potential systematics 
(Galactic extinction, stellar density, synchrotron emis- 
sion) that trace the structure of the Galaxy. Further, 
we trust correlations between the quasars and the LRGs, 
as the LRG sample has already proven to be robust to 
systematic fluctuations. Thus, while we do not consider 
the quasar ACF as a reliable probe of PNG, we will ex- 



ploit the external correlations between the quasars and 
the other data sets. Also in this case, this includes the 
cross-correlation with the CMB, which for the same rea- 
sons should be relatively free from contamination, as also 
confirmed by its fequency independence shown in G12. 



IV. MODELING THE DATA 

A. Data Considered 

We have discussed six different large-scale structure 
data sets, which yield six auto-correlations, fifteen cross- 
correlations and six correlations with the WMAP CMB 
temperature. Our final data set is shown in Fig. |11[ 
including the galaxy-CMB cross-correlations and the 
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Parameter 


Description 


Prior range 


Cosmology 


LJb = ^bh^ 


baryon energy density 


[0.015,0.035] 




ujc = 0.ch 


dark matter energy density 


[0.05,0.2] 




e 


sound horizon at the last-scattering surface 


[0.8,1.2] 




T 


optical depth 


[0.01,0.25] 




log (10^° A,) 


amplitude of primordial perturbations at fcpiv = 0.002/i/Mpc 


[2.9,3.5] 


Our focus 


/nl 


amplitude of skewness-type primordial non-Gaussianity 


[-200, 200] 




9NL 


amplitude of kurtosis-type primordial non-Gaussianity 


[-5- 10^5■ 10"] 




INL 


scale-dependence of bias 


[0,4] 




W 


dark energy equation of state 


[-2.5,-1/3] 


Nuisance 


Asz 


amplitude of the SZ CMB template 


[0,2] 




bh 


6 bias nuisance parameters for each catalog 


[0,3] 




7i 


slope of bias evolution for each catalog 


[0,3] 




Ki 


3 stellar contamination nuisance parameters for each SDSS catalog 


[0,0.1] 






15 dN/dz nuisance parameters for each pair of catalogs 


Gaussian priors 




Qheao 


1 PSF nuisance parameters for the HEAO catalog 


[0.0001,0.001] 


Derived 




dark energy energy density 


— 




To 


age of the Universe 








total matter energy density 






0"8 


LSS amplitude of density fluctuations 






■2re 


redshift of reionization 






Ho 


Hubble constant 





TABLE I. Summary of the parameters used in the likelihood runs, and their prior ranges. We checked that our posterior results 
are never affected by the edges of the priors. 



galaxy-galaxy correlations between all the catalogs. In all 
cases below, we discard the first five angular bins of the 
2MASS-CMB CCF, which is believed to be contaminated 
by the SZ effect. We also discard the zero-lag data points 
in the density-density correlations based on contamina- 
tion from non-linear structure formation, which we are 
not modeling, and from possibly inaccurate shot-noise 
subtraction in the correlations between different catalogs, 
due to the unknown fraction of common sources between 
the catalogs. 

Given the potential impact of systematic contami- 
nants, it is useful to define the following three sub-sets of 
data (we exclude in all cases the 0-lag bins of the density- 
density correlations and the first five data points in the 
2MASS ISW measurement as described above): 

• The 'naive' data set is the full combination, which 
includes all 27 correlation fmictions. 

• The 'fair' sample contains 25 correlations, exclud- 
ing the auto-correlations of NVSS and QSO, which 
are expected to be overall the least reliable parts 
of our data set, due to the high risk or residual 
systematics (as discussed above). 

• The 'conservative' set includes only what we view 
as the most reliable correlations. We exclude all 
the auto-correlations except for the LRG case, as 
auto-correlations are expected to be more prone to 
systematics than the cross-correlations. Further, 
we exclude any data from 2MASS and from the 
main SDSS galaxies; as they are not expected to 
contribute significantly to the PNG measurements 
due to their low bias, we feel it is better to exclude 



them rather than risk any systematic contamina- 
tion that they might bring. Finally, we also exclude 
the LRG-quasar correlation: as both samples are 
derived from SDSS, there is a limited risk of resid- 
ual correlated systematics, particularly related to 
stellar contamination. The resulting set is consists 
of ten correlations only. 

We will discuss the results from these samples, as well 
as the effect of adding the NVSS auto-correlation to the 
'fair' sample, in the following section. Finally, it is worth 
noting that we discard samples by effectively setting the 
covariance for these data to a high value. The likelihood 
estimation from the inverse of the covariance matrix will 
then simply ignore the discarded points. 

In addition to our data, we will consider also the CMB 
data from WMAP7 [3 and in some cases the luminosity 
distances to the Union2 compilation of Type la Super- 
novae |106| . 

We calculate the theoretical correlation functions us- 
ing a modified version of Camb [107] (validated with 
a fully independent code for accuracy). The covariance 
matrix is generated using the Monte Carlo method de- 
scribed in G12, now extended to the whole data set (we 
have n = 351 data points, so this is the dimension of the 
full covariance). We base our Monte Carlo realizations 
on a fiducial flat ACDM model without PNG, which is 
appropriate since we do not find posterior evidence for 
PNG; if such evidence were found, we would need to re- 
calculate our likelihoods with a new covariance matrix, 
based on a non-Gaussian model. To calculate the co- 
variance, we use the linear bias parameters fitted from 
the auto-correlation functions and assumed constant in 
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redshift. 



B. Gaussian Bias and Other Nuisance Parameters 

In addition to the standard cosmological parameters, 
we must include other nuisance parameters to account 
for the systematics and other uncertainties, over which 
we marginaUze. These describe our assumptions about 
the intrinsic bias (in the absence of PNG effects) , uncer- 
tainties in the redshift distributions of the surveys (and 
particularly their overlaps) and other effects. 

For most catalogs i, we assume an intrinsic (Gaussian) 
bias evolution of the form: 



b\iz) 



1 



1 



D-f'{z) 



(21) 



we assume 7^ = 2 in most cases |108j : this bias ap- 
proaches unity in the future, and is therefore appropriate 
for flux-limited samples. However, for the quasars we in- 
stead assume 



6^="(z) 



(22) 



which instead approaches zero in the future. This is 
consistent with observations of quasars }109| IllOj , that 
suggest that the preferred halo mass of quasars is ~ 
10^^ Mq, independent of redshift. Setting 7qso = 1-6 
approximately satisfies this constant-mass condition for 
0.43 < z < 2.4, which spans the median redshifts of 
Refs. |1091 IllOj . and is the span of redshifts at which we 
may reasonably expect to obtain meaningful information 
from the DR6 quasar sample. We are relatively uncer- 
tain of the bias evolution of the NVSS sample, which 
has a broad redshift distribution but also the potential 
to provide strong PNG constraints due to its large sky 
coverage. We therefore have tested allowing a free slope 
7NVSS for this sample and found that this results in neg- 
ligible changes in our constraints. This suggests that our 
results are not strongly affected by uncertainties in the 
bias evolution relationship. 

In addition we allow the possibility to add 15 extra 
parameters /3ij to account for uncertainties in the red- 
shift distributions ipi{z). A factor arises for each pair of 
catalogs, and they are introduced in the density-density 
cross-correlations so that: 



(23) 



Effectively, the uncertainty in ipi{z) becomes more im- 
portant when the overlap between two surveys is smaller. 
For this reason, we introduce a Gaussian prior, with mean 
— {Pij) — 1 and a variance ct| which depends on the 
fiducial survey overlap, such that the r.m.s. is ct^ = 
when the overlap is total (the ACFs), and ap = 2 when 
the overlap is expected to be small. We linearly interpo- 
late the prior between these two values. We also checked 
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FIG. 12. Marginalized posterior probability distributions on a 
choice of cosmological parameters for different models. First 
row: The flat ACDM case. The results are in agreement with 
the concordance model, but we can see that our data (the 
'fair' data set is used) shift the peak of the distribution to- 
wards models with more dark energy. Second row: Same, with 
the addition of local /nl. (The CMB temperature power spec- 
trum can not constrain PNG and is ignored.) The standard 
cosmological parameters do not change, and /nl is consistent 
with zero. Third row: Same, for the case of a wCDM model. 
Here the CMB can not constrain these parameters simulta- 
neously, hence the broad posteriors. The degeneracy between 
flm — M is broken by adding either Type la SNe or our data. 
The HST Ho prior was added to the w run using our data. 



that the results on PNG do not change significantly when 
doubling the r.m.s. of these priors. 

We also add a free stellar contamination fraction Hi for 
the three samples derived from SDSS, used as described 
for the quasars in Eq. 



(20) 



We finally add one extra pa- 
rameter for the HEAO catalog: an amplitude describing 
the clustering-independent auto-correlation due to source 
shot noise spread over the large PSF of the instrument 
aHEAO 



V. RESULTS 

We examine a range of models; starting from a simple 
flat ACDM model, we add primordial non-Gaussianity of 
different types, or the dark energy equation of state w. 
The parameters which we use and their prior ranges are 
shown in Table IH 

We study the cosmological consequences of these data 
via a Monte Carlo Markov Chain (MCMC) method us- 



ing the publicly available Cosmomc code |112j and its 
extension to nested sampling Multinest |113j , which re- 
places the Metropolis-Hastings with a nested sampling al- 
gorithm, achieving in the process higher accuracy, greater 
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FIG. 13. Marginalized posterior probability distributions on a choice of nuisance parameters for different models. First row: 
The bias parameters &o for each catalog. Second row: Same, for the stellar contamination parameters k^, which are introduced 
for the SDSS data only. The wCDM case includes the HST Hq prior. 




Catalog 


A 


± 0"A 


S/N 


expected 


2MASS cut 


1.32 


± 


2.02 


0.7 


0.5 


SDSS gal DR8 


1.29 


± 


0.59 


2.2 


1.6 


SDSS LRG MegaZ 


2.10 


± 


0.84 


2.5 


1.2 


SDSS LRG CMASS 


1.10 


± 


0.61 


1.8 


1.6 


NVSS G12 


1.21 


± 


0.43 


2.8 


2.6 


NVSS 'shuffled' 


1.47 


± 


0.41 


3.6 


2.6 


HEAO 


1.37 


± 


0.57 


2.4 


2.0 


SDSS QSO DR6 


1.46 


± 


0.60 


2.4 


1.7 


TOTAL G12 


1.38 


± 


0.32 


4.4 ± 0.4 


~3.1 


TOTAL this work 


1.39 


± 


0.31 


4.5 ± 0.4 


~3.1 



TABLE II. Single-catalog and total ISW significances when 
updating the LRGs and cleaning the NVSS data. A ~ 0.4 
error on the full signal-to-noise ratio is added as in G12 to 
account for the typical changes resulting from different choices 
in the analysis. 



speed, and an estimation of the Bayesian evidence for 
each model. We ensured to always use high enough accu- 
racy within MuLTiNEST and Camb to obtain numerical 
stability, also enforcing a sufficiently broad k range in the 
integrations. 



A. ACDM 

ISW Significance We first consider a ACDM tem- 
plate with a single amplitude of the ISW A\ as described 
in G08 and G12 this can be done by assuming the CCFs 
in each bin i are written as wj^ = Agi, where gi is the 
template prediction from the fiducial model. Then the 
best-fit A can be calculated from the covariance matrix 
C and from the observed CCFs wj^ by analytically max- 
imizing the likelihood: 

A= ^r ' . (24) 



where in our case p = 13 if considering a single catalog, 
and p = 13 X 6 = 78 for the full combined analysis. The 
variance is given by 

. (25) 

In this way we recover the significances reported in Ta- 
ble [nj which are consistent with our earlier work, and 
with the expected signal-to-noise ratio within the statis- 
tical errors, calculated according to Eq. (6) in G12. The 
significance of the ISW detection in the CMASS LRGs is 
at l.Scr, down from 2.5cr which was found with the MegaZ 
data, as the measured CCF has decreased, in agreement 
with the results by Ref. [S3]; the significance of NVSS af- 
ter the 'shuffling' procedure has increased to 3.6(7, as the 
measured CCF has increased. The overall significance 
has remained nearly unchanged from G12 at 4.5cr. 

Full MCMC We then fully test a simple flat 
ACDM model using Multinest; we can see in the first 
row of Fig. [12] the resulting ID marginalized posteriors 
on a choice of cosmological parameters. Adding our 'fair' 
data to the WMAP7 CMB power spectrum yields results 
which agree within the errors, even though our data pre- 
fer a lower matter density we find 0.21 < fi^ < 0.28 
(95% c.l.) when using our data and the CMB. We show 
in Fig. [13] the marginalized constraints on the most im- 
portant nuisance parameters. We can see that about 5% 
stellar contamination seems present in the main galaxy 
sample, which is reduced to 1% in the quasars, while 
no evidence for residual contamination appears for the 
LRGs, in agreement with our positive assessment of the 
purity of this sample. 

B. Primordial Non-Gaussianity 

The /isTL Model We next consider the local /nl 
model, which contains one extra parameter. When es- 



16 



0.3 



0.25- 



0.2- 



Conservative 
\ 


\ Naive 









-40 -20 20 40 60 80 



0.85 




0.75 



FIG. 14. Marginalized posterior probability distributions for 
/nl versus QrujCrs for different data sets. The naive result 
obtained using the complete data set (red contours, at 68 
and 95% c.L), which would suggest the presence of signifi- 
cant PNG, is not stable. When using only the most reliable 
parts of our compilation, we obtain the conservative result 
(green), which remains stable when adding back most of the 
data (blue), except the quasar and NVSS ACFs. As these 
are the least reliable data, as we discussed in our systematic 
section, we decide to discard them. 



timating the posterior likelihood of /nl from the whole 
'naive' data set, including all auto-correlations, we find 
a strong positive detection: 30 < /nl < 62 at 95% c.L, 
as shown in 2D by the red contours of Fig. [14] This re- 
sult is consistent with previous measurements that used 
similar data sets '36|, but our uncertainty on /nl is re- 
duced due to the fact that we: (i) Effectively include all 
multipoles / (for SDSS data, we include the k parameters 
to marginalize over the known low-/ issues); (ii) Include 
fainter NVSS data (Ref. [36^ applied a 10 mJy cut; we 
apply shuffling to take care of the systematics instead); 
(iii) Include 33% more area for the LRGs; (iv) Include 
additional data sets. 

To assess whether this result is robust, we study how 
it is changed when less reliable data are excluded. When 



Data 


/nl 95% interval 


Conservative 
Fair 

Fair -t- NVSS ACF 
Naive 


-37 < /nl < 25 
-29 < /nl < 31 
-28 < /nl < 57 
30 < /nl < 62 


conservative, witn Pij, tr^ — mz 
Conservative; with /5ij, a™'"' = ln4 


— OO < JNL < ^( 

-37 < /nl < 25 


LRG-LRG 
NVSS-NVSS 
QSO-QSO 


-90 < /nl < 120 
-68 < /nl < 113 
-15 < /nl < 98 



TABLE III. Summary of the /nl constraints. In the top sec- 
tion we show the results from the combined analyses. The 
two runs with the nuisance parameters Pij have two different 
choices for their Gaussian priors. The runs with a single ACF 
below include marginalizations over one bias parameter and 
(for the SDSS catalogs) one stellar contamination parameter. 
The reason the result with the NVSS ACF appears weaker 
than expected is that it features a double peak in /nl. 



we repeat this test with the 'conservative' data, we find 
no detection: —37 < /nl < 25 at 95% c.L, corresponding 



to the green contours of Fig. 14 When we examine the 
'fair' data, we still see no hint for PNG (blue contours 
in Fig. fl4]) 
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we 



By looking at both panels of Fig 
can see no evidence of strong degeneracies between /nl 
and f2m or erg. The results are in agreement with the 
simplest Gaussian ACDM model: the marginalized 95% 
c.L interval from our 'fair' sample is —29 < /nl < 31. 

We found that including also the NVSS AGF leads to a 
higher inferred /nl because of its high AGF in the angu- 
lar range 2 < < 5 deg. In this case, although we include 
more data than in the conservative case, we actually see 
an increase in the size of the errors: —28 < /nl < 57 
(95%). We have found the reason by looking at the re- 
sult when using the NVSS AGF only. The marginal- 
ized likelihood of /nl is in this case broader and, while 
peaked around /nl — 50, it has a long tail extending un- 
til /nl — —50 with a small secondary peak. This means 
that the previous best-fit region around /nl — does not 
get discarded. This suggests that, even after our clean- 
ing procedure, the clustering present in the NVSS AGF 
is still in tension with the other data sets. 

We can see the corresponding marginalized ID distri- 
butions in the second row of Fig. |12| by comparing with 
the first row, we can also notice that best-fit regions of a 
selection of other cosmological parameters do not change 
when adding /nl^ this indicates overall low degeneracy 
between these parameters and /nl- 

When considering the 'naive' data set, we examined 
the importance of the degeneracies between /nl and 
other parameters, particularly the stellar contamination 
K and bias 69. We found that, perhaps surprisingly, there 
is no significant degeneracy with the quasars' k, while the 
degeneracy with 69 is more pronounced: by raising it to 
&o — 1, values as low as /nl — 30 are allowed. When 
using the quasars' AGF however, the Gaussian limit is 
always excluded at > 2a. We found that a stronger de- 
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FIG. 15. Comparison of the marginalized posterior probabil- 
ity distribution on /nl using the parts of our data set giving 
the strongest contributions. We show the results from sin- 
gle cross-correlation functions (top, green), auto-correlations 
(center, blue), and from combined sub-samples of the whole 
data set (bottom, red). The lines correspond to 68 and 95% 
ranges, have been marginalized over the cosmological param- 
eters, and include the WMAP7 CMB priors. The points 
represent the mean values of the posterior likelihoods. The 
results from single auto-correlation functions have also been 
marginalized over one bias parameter and one stellar contam- 
ination fraction (for the SDSS samples). The NVSS ACF 
result appears weaker than expected beacuse it features a 
double peak in /nl. To best present the relative constrain- 
ing power of the cross-correlation measurements, we have 
placed priors on the bias and stellar contamination parame- 
ters, which significantly overstate the constraints these cross- 
correlation allow on their own. See the main text for more 
details. 



generacy between k and /nl is present only when using 
the quasar ACF alone. 

We summarize the constraints on /nl in Table [TTI| and 
in Fig. [15] for clarity. Here we compare the marginal- 
ized results obtained when using the most constrain- 
ing parts of our data set. We can see once again 
that most resuhs agree with Gaussian initial conditions, 
and with each other. When considering single auto- 
correlation functions, we marginalize over cosmology in- 
cluding the WMAP CMB likelihood, and over one bias 
parameter and one stellar contamination fraction (for the 
data derived from SDSS). To better interpret the cross- 
correlations on their own, we have assigned Gaussian pri- 
ors on the relevant bias and stellar contamination pa- 
rameters equal to the posteriors on these parameters ob- 
tained from the 'fair' data. Applying these priors allows 



us to accurately portray the relative importance of each 
cross-correlation to our bottom-line results. Further- 
more, we found that applying the bias prior to the auto- 
correlations would increase the precision of their /nl con- 
straints by a factor of two. Accounting for this factor, the 
LRG auto-correlation is the best-constrained measure- 
ment that enters the 'conservative' data set. When using 
the LRG ACF only we recover a result consistent with the 
recent analysis by Ref. [37], who found —45 < /nl < 195 
at 95% using the spectroscopic sample of the CMASS 
LRGs, which contains ~ 1/3 of the photoz sample we 
use. 

Notice that the factor (61 — 1) within the bias correction 
Ab is the leading contribution that determines the size 
of the /nl error bars. For this reason, the low-bias data 
from 2MASS, the SDSS main galaxies, and HEAO bring 
little information on /nl- Also the external correlations 
of the quasars bring less contribution than it may be 
expected, since the quasar bias at low redshift is also 
low. This explains why the strongest constraints come 
from NVSS, the LRGs and their external correlations. 
For this reason, we have also checked the effect of the 
assumed NVSS bias evolution with one additional run 
where the evolution parameter 7nvss is let free, and we 
found no significant changes in the results. 

The aisjL Model We then extend our model to gen- 
eralized PNG defined in Eq. (|8]): in addition to /nl, 
we thus allow for scale dependence of the bias of any 
slope cnl, which reduces to onl = 2 in the local, scale- 
independent case. We show our marginalized posterior 
likelihood distribution in the top panel of Fig. [THj where 
we can see that, in line with the lack of evidence for /nl, 
there is no evidence for onl either. The full marginalized 
upper limit we find is onl < 1-7 at 95%, but it must be 
born in mind that there is an infinite degeneracy along 
the direction /nl = by construction: thus, this result 
is strongly dependent on our adopted priors, rather than 
being a "stand-alone measurement" . The correspondent 
bound on Uf^^^ can be found using Eq. (10). 

The gNL Model We finally consider the (/nl model. 
We shall here make the optimistic assumption that the 
fitting formula of Eq. ([7| is a reasonable approximation 
to the effect of ^nl , keeping in mind that this may not be 
accurate in all cases due to the low bias of our catalogs. 
Under this assumption we find —4.5-10^ < (?nl < 1.6-10^ 
(95%) if assuming /nl = 0. However as shown by 
Refs. [m ES], and as clear from Eq. (ij), there is a de- 
generacy between /nl and ^nl, as both parameters pro- 
duce a scale dependence of the bias of the same order 
^ fc~^; the degeneracy is alleviated by the different red- 
shift dependences. This is indeed what happens when 
we consider the complete model where both parameters 
are left free: we can see in the bottom panel of Fig. [16] 
that the marginalized posterior presents this degeneracy, 
as demonstrated with A^-body simulations by Ref. |65j . 
Also in this case the Gaussian model remains well within 
the 95% region: the marginalized constraints on the two 
parameters are marginally degraded to —23 < /nl < 42 



18 



2.5 




FIG. 16. Marginalized posterior probability distributions for 
extended PNG models. Top panel: The parameter hnl de- 
fines the scale-dependence of the bias and it is onl = 2 in 
the local, scale-independent model. The two shaded contours 
represent the 95 and 99% confidence regions. There is a ver- 
tical infinite degeneracy along the /nl = 2 direction, which 
is only partially visible due to finite sampling resolution and 
smoothing. Bottom panel: Marginalized constraints on the 
/nl — 5nl degeneracy (68 and 95% regions). As both pa- 
rameters produce the same scale dependence of the bias, they 
are degenerate, but only partially, as the redshift dependence 
is different. Note that the (Jnl constraints are optimistic, as 
they assume the validity of the fitting formula by Ref. for 
our data. 



and -5.8-10^ < 5nl < 1.7-10^ at 95% respectively when 
they are both set free. 



C. Dark Energy 
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FIG. 17. Marginalized posterior probability distributions for 
rim vs w for difi'erent data sets: CMB only (blue), CMB -f 
Type la Supernovae (red), and CMB -I- our data (green). The 
results are in agreement with the concordance ACDM model. 
The EST Hq prior was added to the run using our data. 



which is broken instead by either Type la Supernovae 
or our data. This can be seen more clearly in the 2D 
plot of Fig. [iTj Here we show that the results from our 
data and the CMB power spectrum reduce the parameter 
space significantly around the concordance values. Notice 
that we have added the HST prior |114j on the Hubble 
constant Hq — 74.2 ± 3.6, since we found that other- 
wise significant degeneracies appear between w, the bias 
parameters, and iJo, due to compensating effects in the 
growth function D{z). When marginalizing over all other 
parameters, we find from our data —1.07 < w < —0.69 
at 95%. This confirms the overall agreement with the 
standard ACDM paradigm of the ISW measurements. 



Without the ISW effect from the presence of dark 
energy, we would not be able to use the CMB cross- 
correlations to constrain PNG. In the above analysis, we 
jointly fitted the dark energy density and PNG to test the 
sensitivity of the PNG measurements to our assumptions 
about DE. 

The wCDM Model We next also consider models 
with Gaussian initial conditions, and use the ISW data 
to constrain the simplest dynamical dark energy of equa- 
tion of state w. We present in the third row of Fig. [12] 
the marginalized ID posterior probability distributions 
when considering this wCDM model. Here we can see 
that, as it is well known, the CMB temperature power 
spectrum alone can not break the fi^ — w degeneracy, 



D. Model Selection 

Nested sampling allows model selection in addition to 
parameter estimation, thanks to the calculation of the 
Bayesian evidence factor for each model. Briefly, if we 
assume a model M of parameters 0, and we compare it 
with data D, Bayes' theorem states that the posterior 
probability distribution V on the parameters is given by 



^(e)H(9) 
Z{M) 



(26) 



where the prior is H(0) — P(8|M), the likelihood 
is C{Q) — P(D|8,Af), and the Bayesian evidence is 
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Parameters 


Data 


In(^) 


InS = Aln(2) 


Odds 


Interpretation 


ACDM 


CMB + our 


-3983.41 ±0.13 





— 


— 


+ /nl 


('fair') 


-3985.03 ±0.13 


-1.62 ±0.18 


1 : 5 


Weak evidence against /nl 


+ 5NL 




-3985.31 ±0.13 


-1.90 ±0.18 


1 : 7 


Weak evidence against (/nl 


+ /nl + 5NL 




-3986.59 ±0.13 


-3.18 ±0.19 


1 : 24 


Moderate evidence against /nl + 3nl 


wCDM 




—3984.31 ± 0.13 


—0.90 ± 0.18 


2 : 5 


Inconclusive 


ACDM 


CMB alone 


-3751.37 ±0.08 









wCDM 




-3751.84 ±0.09 


-0.47 ±0.12 


5 : 8 


Inconclusive 


ACDM 


CMB + SN 


-4017.15 ±0.09 









wCDM 




-4019.10 ±0.09 


-1.95 ±0.13 


1 : 7 


Weak evidence against w 



TABLE IV. Bayesian model selection. For each data set we compare models with extra parameters with the baseline ACDM. 
We can see that our 'fair' data in combination with the CMB provide weak evidence against models with one extra PNG 
parameter, and moderate evidence against a more complex model with both /nl,5nl present. 



Z{M) = P{D\M). When interested in parameter esti- 
mation only, it is common to neglect the evidence Z and 
simply study the posteriors of arbitrary normalization for 
a given model. However here we will also compare differ- 
ent models (e.g. with and without PNG), and for this we 
will use the evidence. Model selection between two mod- 
els M, can be performed by calculating the ratio of 
their evidences, also called Bayes factor: B = Z^/Zm, 
so that InB = AlnZ. This factor, which incorporates 
Occam's razor by penalizing models with unnecessary ex- 
tra parameters, expresses the odds between the models 
given the data, and can be qualitatively interpreted with 
the heuristic Jeffrey's scale, which states that the model 
selection is inconclusive if | In i? | < 1 , and that the ev- 
idence for one model is weak, moderate, or strong for 
InB > 1, InS > 2.5, \nB > 5 respectively pT5] . 

To test non-Gaussianity, we can compare the Bayesian 
evidence Z of the models with and without PNG. As 
shown in Table|TVj the Bayes factor between these models 
is AlnZ = —1.62 in the case of /nl and AlnZ = —1.90 
in the case of (7nl , which are both interpreted as weak ev- 
idence against these models according to Jeffrey's scale, 
as the two models have odds ~ 1 : 5 and 1 : 7 compared 
with ACDM respectively given our data [115 ! . The situ- 
ation is more defined for the extended model with both 
/nl and ^NL^ in this case we find AlnZ = —3.18, cor- 
responding to odds of 1 : 24, which is seen as moderate 
evidence against this model. 

We also performed a Bayesian model selection for dark 
energy: as shown in Table |IV[ when using the CMB and 
our data this test is inconclusive: AlnZ = —0.90, which 
means the odds are 2 : 5 when comparing ACDM with 
a wCDM model. The result becomes more defined only 
when combining CMB and data from Type la Super- 
novae: in this case AlnZ = —1.95, interpreted as weak 
evidence against the wCDM model. 

In general it has to be kept in mind that the Bayesian 
selection is affected by the assumed choice of priors, as 
the prior normalization implies the evidence is decreased 
if the prior range is broader. We have however chosen a 
relatively narrow prior on /nl and ^nl compared with 
existing constraints, meaning we are not unnecessarily 
penalizing the evidence for these PNG models. 



VI. CONCLUSIONS 

Scale-dependent bias has arisen as a key means of de- 
tecting primordial non-Gaussianity, but its detection can 
be severely compromised by a number of potential sys- 
tematics. These, if not accounted for, can introduce 
spurious large-scale structure and bias the inferred con- 
straints on PNG. We have therefore conducted an analy- 
sis that carefully accounts for these systematic concerns 
and maximizes the information obtained from the most 
robust sources, in order to obtain the most precise to- 
date PNG measurements. 

In particular, large differences are observed when com- 
paring the MegaZ and CMASS LRG samples. As shown 
in Ref. |53j , these appear to arise mainly from systematic 
relationships with foreground stars which, if unaccounted 
for, cause spurious large-scale power. Similarly, correct- 
ing the known declination-dependent systematic in NVSS 
significantly reduces the inferred large-scale structure. 
Even so, some large-scale power remains in NVSS, but 
we have shown that this is potentially due to a simi- 
lar systematic connected to variations in the mean den- 
sity as a function of the right ascension. Finally, while 
quasar samples are potentially ideal for searching for pri- 
mordial non-Gaussianity, we have confirmed earlier work 
[5 9) showing that the quasar auto-correlation measure- 
ments are compromised by systematic errors. In particu- 
lar, there are correlations between the mean quasar den- 
sity and a number of systematics including stellar density, 
galactic reddening and sky brightness. 

While the existence of such systematics can be discour- 
aging, they can often be corrected when they are under- 
stood. In addition, cross-correlations between different 
data sets, which should generally contain different sys- 
tematics, offer a much more robust means of addressing 
the question of primordial non-Gaussianity. 

One example of this can be seen in the analysis of cross- 
correlations with the CMB, previously analysed in G12. 
Comparing the MegaZ and CMASS measurements, the 
CMB cross-correlation was seen to be much less affected 
than the auto-correlations. While there were differences, 
these are partially accounted for by the larger sky cov- 
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erage of CMASS; the extra area has associated cosmic 
variance and is potentiaUy more subject to extinction 
systematics. The signal decreased with respect to the 
G12 resuh, but the error bars are also smaller due to the 
larger DR8 sky coverage: the total signal-to-noise of this 
ISW detection alone is at the 1.8cr level, in agreement 
with other recent measurements [54] . Similarly, account- 
ing for the newly discovered NVSS systematic changes 
its cross-correlation somewhat but not dramatically, ac- 
tually raising the inferred signal by 20%. However, the 
total ISW significance from the combined data is in the 
end unchanged with respect to G12, as the increase in 
NVSS compensates the decrease in the LRGs. 

While measurements of cross-correlations are more ro- 
bust, they require more information for their interpreta- 
tion, in particular understanding how the surveys overlap 
in redshift. We treat these uncertainties by marginalizing 
over nuisance parameters associated with the strength 
of the cross-corrrelation, so that it is its shape, rather 
than its amplitude, that constrains the PNG. Effec- 
tively, the overlap is determined by the small scale cross- 
correlations, allowing us to use the large-scale measure- 
ments to constrain PNG. Fundamentally, it is the lack of 
large-scale excess power in any of the cross-correlations 
which is responsible for our strong constraints. The pre- 
cise strength of the constraints will depend on our as- 
sumptions about the bias evolution, particularly for the 
deepest surveys; however, the bias models which we use 
are consistent with our physical understanding of the 
tracers and other observations. Further, we have applied 
a general bias model, one that effectively assumes only 
that the bias smoothly tends towards unity as the redshift 
decreases, to the NVSS sample and found no significant 
change in our results. 

We have measured that the allowed amount of PNG 
is —37 < /nl < 25 at 95%. This bound was determined 
from our conservative analysis, for which we excluded 
all auto-correlations apart from the CMASS LRGs, and 
instead relied on the the cross-correlations between large- 
scale structure data sets and also CMB cross-correlations 
which are sensitive to scale-dependent bias through the 
ISW effect. We found the data consistent with the sim- 
plest ACDM model in all cases. This is significantly 
different than the positive signal detected if all possi- 
ble auto-correlations are included, as the signal is driven 
up by the large-scale clustering in the quasar catalog. 
Including the NVSS auto-correlation actually broadens 
the error bar compared to the conservative result, as the 
survey has some residual large-scale power even after the 
corrections for declination and right ascension dependent 
mean density. 

Our measurements have also allowed us to constrain 
more general PNG models, as well as models of dark 
energy. As yet, there are no indications of deviations 
from Gaussianity, or dynamic dark energy. Imminent re- 
sults from Planck and future large-scale structure mea- 
surements, such as will arise from DES, LSST and Eu- 
clid, may well determine whether the simplest infiation- 



ary models are correct. However, the more sensitive the 
results, the more prone they are to potential systematics, 
making a careful comparison and correlation between a 
range of different measurements indispensable. 
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Appendix A: Fast Calculation of PNG Spectra 

Our analysis requires scanning over a large number 
of possible model parameters, which is potentially very 
computationally expensive if a CMB code such as Camb 
[107] must be rerun for every parameter choice. Here we 
describe how we can make the calculation of the PNG 
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effects more efBcient, so that a range of bias choices and 
/nl values can be calculated using a single run of Camb. 



1. Gaussian Case 

In the Gaussian case, number density perturbations at 
a position x and redshift z are defined at linear level as 



^g(x,2;) = 6g(z)5(x,z), 



(Al) 



where hg{z) is the galaxy bias and <5(x,z) are the under- 
lying dark matter perturbations. Projection onto the ce- 
lestial sphere is achieved by summing over the past light 
cone, so that the 2D number density perturbation in a 
direction fi is 

5g(n) = J dz(p{z)b{z)S{:>i,0)D{z) , (A2) 

where D{z) is the linear growth function normalized to 
unity today and ip{z) is the normalized redshift distribu- 
tion of the sources. 

We can expand as usual these perturbations into har- 
monic space by means of a Fourier transform and an ex- 
pansion in spherical harmonics Yim 



(A3) 



k /,m 

where ji are the spherical Bessel functions. We thus get 
the harmonic coefRcients of the perturbations 



ai 



= ^W^(k,z)y4(r!fc) / dz^{z)bg{z)D{zMkr). 
k •' 

We then introduce for simplicity the quantity 

fi{k) = [ dz^{z)b{z)D{z)ji{kr), (A4) 



so finally the projected angular power spectrum of the 
galaxy density perturbations can be written in a compact 
form in function of these quantities as 



Cf« ^ (a,,„aL) ^ - I dkP P{k) ff{k) . (A5) 



2. Non-Gaussian Case 

When local non-Gaussianity of the /nl type is intro- 
duced, we know that at linear level the galactic bias will 
be corrected by 



M{k,z) 



2/NL^c[bl(^)-l] 

a{k, z) 



(A6) 



where we have used the definition of a of Eq. ([s]) . While 
redshift evolution of the bias is difficult to model, a gen- 
eral template can be written as a function of three pa- 
rameters &0 7^-i,7 |108j as 



Df{z) 



(A7) 



We will here consider three simplified cases, introducing 
different templates for it, such as simply a constant 



bi{z) = bo , 



hi{z) = 1-1- 



D-^{zy 



(A8) 



(A9) 



which satisfies the property of approaching unity in the 
future, and is therefore appropriate for flux-limited sam- 
ples, and 



D-^iz) 



(AlO) 



which instead tends to zero in the future, consistent with 
observations of quasars |109i 1110) . 



We can now generalize the 2D spectrum of Eq. ( A5 ) 



between two galaxy catalogs i, j by noticing how the bias 
alterations due to PNG can be simply absorbed into the 
source function which can be split in two: 



dk k^P{k) 



(All) 



flik) + 



9l{k) 



a{k,0) 



ink) 



a{k,0) 



where we have absorbed the redshift dependence of a (a 
factor of D) into the integrand of , and 



//(fc) = j dz^,{z)b\{z)D{z)ji{kr) (A12) 
g\{k) = fdz ^,{z) {2 /nl 5c [h\{z) - l] } 3i{kr) . 



In the integrand of g] we have simplified the factor of 
D{z) from a, so that its original D{z) has disappeared. 



3. Fast Decomposition of the Spectra 

It is useful to decompose the power spectra in a way 
which allows us to extract the bias and the /nl param- 
eters outside of the integrals, so that the calculation is 
faster. 
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a. Bias Parameterization h\ (2:) = 1 + JyTJJ) 



b. Bias Parameterization bi{z) — 



D-i(z) 



To fix the ideas, let us assume for the redshift evolution 
of the bias that for a catalog i it holds 



— 1 ^* 
61(z)-l= " - ^° 



(A13) 



We want now to express the power spectra using the pa- 
rameter 6^Q. The sources can be written explicitly as 



flik)^ I dz^,{z)D{z)ji{kr) 
-V 



(A14) 



dz „';°'^fl , ji{kr) 



Df-'^{z) ■ 



gl {k)^2 /nl 5, J dz 1^ ji (kr). 
We can introduce here three source terms 



dz (p^{z) D{z)ji{kr) 



9{ 



dz 
2 5r I dz 



Df^iz) 



ji{kr) 
ji{kr) 



(A15) 



which are defined inside Camb in the file equations . f 90, 
subroutine output. The previous expressions can be sim- 
ply rewritten as 



fiik) = ft'ik)+bufi'\k) 

gtik) = hLblo9l'\k), 



(A16) 



which can be done in cmbmain.f90, subroutine 
calcScalCls. In this way the dependencies on both pa- 
rameters and /nl are factored out; the power spectra 
can then be easily calculated using Eq. (All I. 



In this case, b\{z) - 1 = 5^^^ - 1 = jj^i^) - 1' so 
that we can write (notice here we use bo as nuisance pa- 
rameter) : 



flik)=bl, J dz^^^jiikr) 

= -2/nl<5c j dzLpi{z)ji{kr) 

+ 2 bl /nl 5c J dz 1^ ji (kr) . (A17) 



We can define then the source terms 



fi'^ = Jdz J^l^f^^^ Jiikr) 



g\'° = -5c J dz (pi{z) ii{kr) 

g:-'^25cJdz^^Mkr), 
so the previous expressions can be rewritten as 



(A18) 



9m^hL9l''ik) + hLbl,gl'\k). (A19) 



The constant bias case 61(2:) — bo can be also derived 
from these expressions. 
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